import numpy
import matplotlib.pyplot as mpt
import math



#输入为list类型参数，生育率和存活率列表
def init_Leslie(birth,save):
    #Leslie矩阵
    listL = [[0 for k in range(20)] for m in range(20)]
    for i in range(20):
        listL[0][i] = birth[i]
    for i in range(19):
        listL[i+1][i] = save[i]
    L = numpy.matrix(listL)
    return L

#n为待预测年份，（每5年预测一次），输入为Leslie矩阵和初始人口列向量
def run_Leslie(L,save,people,n):
    for i in range(n//5):
        people = numpy.matmul(L,people)
        people[-1] = people[-1] * save[-1]
    return people

#birth各年龄段生育率 5岁唯一档
#save各年龄段存活率（1-死亡率）
# people各个年龄段初始人口
# n 为运行年数，（5的倍数）
def run_module(n,people,birth,save):
    people = numpy.matrix(people,numpy.float32)
    people = people.transpose()
    #运行程序
    L = init_Leslie(birth,save)
    result_people = run_Leslie(L,save,people,n)

    return result_people


#求老年人口比重
def get_sum(result_people,people):
    result_people = result_people.tolist()
    result_people_list=[]
    for i in range(len(result_people)):
        result_people_list.append(result_people[i][0])
    return math.fsum(result_people_list[13:-1])/math.fsum(result_people_list[0:-1])


def main(birth_path,people_path):
    # 计算死亡率，死亡率取2002年，单位为千分之
    save = []
    death = [6.04, 0.55, 0.42, 0.63, 0.97, 1.11, 1.33, 1.68, 2.40, 3.47, 5.49, 8.68, 14.79, 24.37, 42.29, 66.49, 111.68,
             160.30, 242.06, 300.70]
    for i in range(len(death)):
        death[i] = float(death[i])/ 1000
        save.append((1 - death[i])**5)

    # 计算出生率
    birth1 = []
    birth2 = []
    birth3 = []
    with open(birth_path, 'r', encoding='utf-8') as f:
        for line in f.readlines():
            sum = 0
            listline = line.split('\t')
            for i in listline:
                sum += float(i)
            E = sum / (len(line.split('\t')))
            birth1.append(E * 5 / 1000)
            # 出生率修正，即比15年增长多少
            birth2.append(E * 5 * 1.2 / 1000)
            birth3.append(E * 5 * 1.05 / 1000)

    # 人口数
    people_linshi = []
    with open(people_path, 'r', encoding='utf-8') as f:
        for line in f.readlines():
            listline = line.split('\t')
            people_line = []
            for item in listline:
                people_line.append(float(item))
            people_linshi.append(people_line)

    people = [[0 for i in range(len(people_linshi[0]))] for j in range(len(people_linshi) - 1)]
    for i in range(len(people_linshi) - 1):
        for j in range(len(people_linshi[i])):
            people[i][j] = people_linshi[i][j] / people_linshi[-1][j]
    # 获取2001年人口数目
    people_l = []
    for i in range(len(people)):
        people_l.append(people[i][0])
    # 运算,2003年开放二孩，2016年开放二胎
    listyear = [i * 5 + 2001 for i in range(20)]
    listsumpeople = []
    for i in range(20):
        if (i <= 1):
            result_ = run_module(5 * i, people_l, birth1, save)
            listsumpeople.append(get_sum(result_, people))
        elif (1 <i <= 3):
            result_ = run_module(5 * i, people_l, birth3, save)
            listsumpeople.append(get_sum(result_, people))
        else:
            result_ = run_module(5 * i, people_l, birth2, save)
            listsumpeople.append(get_sum(result_, people))
    # 画图
    mpt.xlabel("Year")
    mpt.ylabel("Percentage of People Over 65")
    mpt.plot(numpy.array(listyear), numpy.array(listsumpeople),label="Allow Two Child Policy in 2006 and Twice Birth in 2016")
    mpt.legend(loc="upper right",ncol=1)
    mpt.title("Data in Current Policy with Two Child Policy in 2006 and Twice Birth in 2016")
    mpt.show()

if __name__ =="__main__":
    main("C:\\Users\\32628\\Desktop\\birth.txt","C:\\Users\\32628\\Desktop\\people.txt")